Rappel contexte

Design, contrastes, et test DE Quelle différence avec une seule condition?

EgdeR and DESeq2 Fit a negative binomial generalized log-linear model to the read counts for each gene. Conduct genewise statistical tests for a given coefficient or coefficient contrast. Ce modèle est de la forme

\(log(\mu_{gj}) \hat= \beta_o + \beta X\) knowming that \(X \sim NB(\mu_{gj}, \sigma_g)\)

Import des données

data <- read.csv("quantifFiles/quantifGenes.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, grepl("R", colnames(data))]
head(data)

            R3   R2   R1   R5  R4   R6  R7  R11  R10  R12  R13  R17  R16  R14
AT1G01010 1526 1006 1116 1275 967 1018 854 1132 1294 1364 2325 2113 2193 2564
AT1G01020  416  285  289  349 364  370 226  513  502  561  461  407  432  614
AT1G01030   31   15   19   29  36   28  12   47   34   47   18   40   32   44
           R15  R18  R24  R19  R22  R23  R21  R20  R9  R8
AT1G01010 2364 2074 1987 2027 1754 1697 1537 1898 816 912
AT1G01020  380  502  484  467  426  415  413  462 223 312
AT1G01030   37   27   42   39   36   40   37   37  15  19
 [ reached 'max' / getOption("max.print") -- omitted 3 rows ]

dim(data)

[1] 28775    24

annot <- read.csv("Code_for_RNAseq_CO2_N_Fr.csv", h = T, sep = ";")
conditions <- as.vector(unique(annot$Sample))
annot$ID <- paste0("R", annot$Code)

annot$condition <- substr(conditions, 1, nchar(conditions) - 1)
cond <- unique(substr(conditions, 1, nchar(conditions) - 1))
cond <- cond[grepl("At", cond)]

getCondition <- function(id) {
    return(annot[annot$ID == id, "condition"])
}

getExactCondition <- function(id) {
    return(annot[annot$ID == id, "Sample"])
}

getLabel <- function(id) {
    text <- as.vector(annot[annot$ID == id, "Sample"])
    res <- ""
    nb <- substr(text, nchar(text), nchar(text))
    if (grepl("Ambient", text)) {
        res = paste0(res, "c")
    } else {
        res = paste0(res, "C")
    }
    if (grepl("High", text)) {
        res = paste0(res, "N")
    } else {
        res = paste0(res, "n")
    }
    if (grepl("Starv", text)) {
        res = paste0(res, "f")
    } else {
        res = paste0(res, "F")
    }
    res = paste0(res, "_", nb)
    return(res)
}

getLabel("R6")

[1] "cnF_3"

getExactCondition("R4")

[1] At_AmbientCO2_LowNitrate_Fe1
48 Levels: At_AmbientCO2_HighNitrate_Fe1 ... Sl_ElevatedCO2_LowNitrate_FeStarvation3

TCC (Tag Count Comparison)

TCC is an R package that provides a series of functions for differential expression analysis of tag count data. The package incorporates multi-step normalization methods, whose strategy is to remove potential DEGs before performing the data normalization. The normalization function based on this DEG elimination strategy (DEGES) includes (i) the original TbT method based on DEGES for two-group data with or without replicates, (ii) much faster methods for two-group data with or without replicates, and (iii) methods for multi-group comparison. TCC provides a simple unified interface to perform such analyses with combinations of functions provided by edgeR, DESeq, and baySeq.

On crée l’objet TCC avec le design souhaité, et on filtre les gènes avec de faibles expressions (paramètre low.count).

Lors de la normalisation (DEGES,iedgeR), on fait un premier calcul des gènes DE, pour pouvoir les enlever lors de la normalisation. Le maramètre test.method permet de choisir la manière de détecter les genes DE (edgeR, DEsqe2, ou tBt (très long)). On peut répéter cette procédure jusqu’à la convergence des facteurs de taille des librairies, d’ou le i.

# design
keep <- rowSums(data) >= 10
data <- data[keep, ]
group <- sapply(colnames(data), getCondition)
colnames(data) <- sapply(colnames(data), getLabel)
tcc <- new("TCC", data, group)
tcc

Count:
          cNF_3 cNF_2 cNF_1 cnF_2 cnF_1 cnF_3 CNF_1 CnF_2 CnF_1 CnF_3 cNf_1
AT1G01010  1526  1006  1116  1275   967  1018   854  1132  1294  1364  2325
AT1G01020   416   285   289   349   364   370   226   513   502   561   461
AT1G01030    31    15    19    29    36    28    12    47    34    47    18
          cnf_2 cnf_1 cNf_2 cNf_3 cnf_3 Cnf_3 CNf_1 Cnf_1 Cnf_2 CNf_3 CNf_2
AT1G01010  2113  2193  2564  2364  2074  1987  2027  1754  1697  1537  1898
AT1G01020   407   432   614   380   502   484   467   426   415   413   462
AT1G01030    40    32    44    37    27    42    39    36    40    37    37
          CNF_3 CNF_2
AT1G01010   816   912
AT1G01020   223   312
AT1G01030    15    19
 [ reached getOption("max.print") -- omitted 3 rows ]

Sample:
                                        group norm.factors lib.sizes
cNF_3            At_AmbientCO2_HighNitrate_Fe            1  32032772
cNF_2            At_AmbientCO2_HighNitrate_Fe            1  25708325
cNF_1            At_AmbientCO2_HighNitrate_Fe            1  25388926
cnF_2             At_AmbientCO2_LowNitrate_Fe            1  29410725
cnF_1             At_AmbientCO2_LowNitrate_Fe            1  26720600
cnF_3             At_AmbientCO2_LowNitrate_Fe            1  26429519
CNF_1           At_ElevatedCO2_HighNitrate_Fe            1  16944203
CnF_2            At_ElevatedCO2_LowNitrate_Fe            1  28518287
CnF_1            At_ElevatedCO2_LowNitrate_Fe            1  30134658
CnF_3            At_ElevatedCO2_LowNitrate_Fe            1  30939161
cNf_1  At_AmbientCO2_HighNitrate_FeStarvation            1  30097179
cnf_2   At_AmbientCO2_LowNitrate_FeStarvation            1  29411029
cnf_1   At_AmbientCO2_LowNitrate_FeStarvation            1  32417299
cNf_2  At_AmbientCO2_HighNitrate_FeStarvation            1  34497190
cNf_3  At_AmbientCO2_HighNitrate_FeStarvation            1  31947050
cnf_3   At_AmbientCO2_LowNitrate_FeStarvation            1  30166724
Cnf_3  At_ElevatedCO2_LowNitrate_FeStarvation            1  33780095
CNf_1 At_ElevatedCO2_HighNitrate_FeStarvation            1  30427501
Cnf_1  At_ElevatedCO2_LowNitrate_FeStarvation            1  29577130
Cnf_2  At_ElevatedCO2_LowNitrate_FeStarvation            1  30405383
CNf_3 At_ElevatedCO2_HighNitrate_FeStarvation            1  25049240
CNf_2 At_ElevatedCO2_HighNitrate_FeStarvation            1  28496520
CNF_3           At_ElevatedCO2_HighNitrate_Fe            1  18823627
CNF_2           At_ElevatedCO2_HighNitrate_Fe            1  20700711

tcc <- filterLowCountGenes(tcc, low.count = 10)

# Normalisation
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", iteration = 1, 
    FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors

    cNF_3     cNF_2     cNF_1     cnF_2     cnF_1     cnF_3     CNF_1     CnF_2 
0.9662816 0.9730845 0.8976118 0.9973121 1.0475825 1.0072335 0.9752235 1.0522596 
    CnF_1     CnF_3     cNf_1     cnf_2     cnf_1     cNf_2     cNf_3     cnf_3 
1.0186883 1.0510472 1.0019289 0.9977532 1.0188031 1.0400636 0.9948079 1.0087956 
    Cnf_3     CNf_1     Cnf_1     Cnf_2     CNf_3     CNf_2     CNF_3     CNF_2 
0.9381444 0.9961312 1.0067383 1.0059904 1.0542699 1.0349814 0.9249584 0.9903091 

tcc$DEGES$execution.time

   user  system elapsed 
 11.340   0.797  12.188 

s <- sample(rownames(tcc$count), size = 200)
heatmap(as.matrix(tcc$count[s, ]), main = "Before normalisation")

normalized.count <- getNormalizedData(tcc)
heatmap(as.matrix(normalized.count[s, ]), main = "After normalisation")

# DEtest
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.001)
result <- getResult(tcc, sort = TRUE)
DEgenes <- subset(result, estimatedDEG == 1)
print(paste(dim(DEgenes)[1], " genes DE"))

[1] "13655  genes DE"

top = DEgenes[1:5000, ]
head(result)

    gene_id a.value m.value p.value q.value rank estimatedDEG
1 AT1G02850      NA      NA       0       0   45            1
2 AT1G03870      NA      NA       0       0   45            1
3 AT1G05660      NA      NA       0       0   45            1
4 AT1G06120      NA      NA       0       0   45            1
5 AT1G09932      NA      NA       0       0   45            1
6 AT1G12030      NA      NA       0       0   45            1

plot(tcc, ylim = c(-11, 11), main = "Nitrate effect en conditions normales")

plot(tcc, group = c("At_AmbientCO2_HighNitrate_Fe", "At_ElevatedCO2_HighNitrate_Fe"), 
    main = "Effet CO2 en conditions normales", ylim = c(-11, 11))

plot(tcc, group = c("At_AmbientCO2_HighNitrate_FeStarvation", "At_ElevatedCO2_HighNitrate_FeStarvation"), 
    main = "Effet CO2 en Fe starvation", ylim = c(-11, 11))

plot(tcc, group = c("At_AmbientCO2_LowNitrate_Fe", "At_ElevatedCO2_LowNitrate_Fe"), 
    main = "Effet CO2 en low Nitrate", ylim = c(-11, 11))

plot(tcc, group = c("At_AmbientCO2_HighNitrate_Fe", "At_AmbientCO2_HighNitrate_FeStarvation"), 
    main = "Effet Fe en conditions normales", ylim = c(-11, 11))

plot(tcc, group = c("At_AmbientCO2_LowNitrate_FeStarvation", "At_AmbientCO2_HighNitrate_FeStarvation"), 
    main = "Effet Nitrate en Fe starvation", ylim = c(-11, 11))

mat <- as.matrix(normalized.count[top$gene_id, ])
heatmap(mat)

# ggplot(melted_res, aes(Var2, Var1, fill= log(value))) + geom_tile() +
# scale_fill_distiller(palette = 'RdPu') + theme(axis.text.x = element_text(size
# = 0.5, angle = 320, hjust = 0, colour = 'grey50')) library(viridis) heatmap <-
# ggplot(melted_res, aes(Var2, Var1, fill= log(value))) + geom_raster() +
# scale_fill_viridis(trans='sqrt') + theme(axis.text.x=element_text(angle=65,
# hjust=1), axis.text.y=element_blank(), axis.ticks.y=element_blank()) heatmap

Fonction pour visualiser un gène étant donné son identifiant, en terme d’expression normalisée :

getExpression <- function(gene, conds = "all") {
    print(conds)
    if (conds == "all") {
        conds = colnames(normalized.count)
    } else {
        conds = grepl(conds[1], colnames(normalized.count)) | grepl(conds[2], colnames(normalized.count))
    }
    df <- normalized.count[gene, conds]
    library(reshape2)
    d <- melt(df)
    d$group = str_split_fixed(rownames(d), "_", 2)[, 1]
    
    p <- ggplot(data = d, aes(x = group, y = value, fill = group)) + geom_dotplot(binaxis = "y", 
        stackdir = "center") + theme(axis.text.x = element_text(angle = 320, hjust = 0, 
        colour = "grey50"), plot.title = element_text(size = 14, face = "bold")) + 
        ggtitle(paste("Normalized expression for ", gene))
    p
}

getExpression("AT1G62510")

[1] "all"

getExpression("AT1G62510", conds = c("cNF", "CNF"))

[1] "cNF" "CNF"

getExpression("AT1G73120")

[1] "all"

Il semble que l’effet du fer soit le plus fort, et celui qui amplifie les autres effets. La technique utilisée ici identifie des DEG globalement, sans séparer lesquels sont dûs à quel effet.

Comme les contrastes ne sont pas donnés, on considère toutes les conditions comme 8 conditions indpendantes.

Analyse globale avec edgeR

TCC fait le test globalement, avec un layout “one way”, qui prend les 8 conditions comme toutes différentes. En spécifiant un design de combinatoire \(2*2*2\), on peut avoir des coefficients testés contre 0 pour chacune des variables et leur interactions.

# data
d <- data[, !grepl("cluster", colnames(data))]


# fixation correcte des contrastes et formule
groups = str_split_fixed(colnames(d), "_", 2)[, 1]


y <- DGEList(counts = d, group = groups)
y$samples

      group lib.size norm.factors
cNF_3   cNF 32032772            1
cNF_2   cNF 25708325            1
cNF_1   cNF 25388926            1
cnF_2   cnF 29410725            1
cnF_1   cnF 26720600            1
cnF_3   cnF 26429519            1
CNF_1   CNF 16944203            1
CnF_2   CnF 28518287            1
CnF_1   CnF 30134658            1
CnF_3   CnF 30939161            1
cNf_1   cNf 30097179            1
cnf_2   cnf 29411029            1
cnf_1   cnf 32417299            1
cNf_2   cNf 34497190            1
cNf_3   cNf 31947050            1
cnf_3   cnf 30166724            1
Cnf_3   Cnf 33780095            1
CNf_1   CNf 30427501            1
Cnf_1   Cnf 29577130            1
Cnf_2   Cnf 30405383            1
CNf_3   CNf 25049240            1
CNf_2   CNf 28496520            1
CNF_3   CNF 18823627            1
CNF_2   CNF 20700711            1

# filtre
keep <- filterByExpr(y)
y <- y[keep, keep.lib.sizes = FALSE]

# normalisation
y <- calcNormFactors(y)
norm <- cpm(y, normalized.lib.sizes = TRUE)
not_norm <- cpm(y, normalized.lib.sizes = FALSE)
# rd_genes = sample(rownames(norm), 500) heatmap(not_norm[rd_genes,], main =
# 'Sans normalisation') heatmap(norm[rd_genes,], main = 'Avec normalisation')


# Definition du design et des bonnes conditions de référence

co2 <- str_split_fixed(groups, "", 3)[, 1]
nitrate <- factor(str_split_fixed(groups, "", 3)[, 2])
nitrate <- relevel(nitrate, "N")

fer <- factor(str_split_fixed(groups, "", 3)[, 3])
fer = relevel(fer, "F")

data.frame(Sample = groups, co2, nitrate, fer)

   Sample co2 nitrate fer
1     cNF   c       N   F
2     cNF   c       N   F
3     cNF   c       N   F
4     cnF   c       n   F
5     cnF   c       n   F
6     cnF   c       n   F
7     CNF   C       N   F
8     CnF   C       n   F
9     CnF   C       n   F
10    CnF   C       n   F
11    cNf   c       N   f
12    cnf   c       n   f
13    cnf   c       n   f
14    cNf   c       N   f
15    cNf   c       N   f
16    cnf   c       n   f
17    Cnf   C       n   f
18    CNf   C       N   f
 [ reached 'max' / getOption("max.print") -- omitted 6 rows ]

design <- model.matrix(~co2 * nitrate * fer)
rownames(design) <- colnames(y)

kable(design)

# design <- model.matrix(~0 + groups) con <- makeContrasts((groupscnF -
# groupscnf), levels=design) fit <- glmQLFit(y, contrast = con) qlf.2vs1 <-
# glmQLFTest(fit) a = topTags(qlf.2vs1, n= 1000) a$table

# estimation of dispersion and tests
y <- edgeR::estimateDisp(y, design)
y$common.dispersion

[1] 0.01113512

plotBCV(y)

fit <- glmFit(y, design)

colnames(design)

[1] "(Intercept)"        "co2C"               "nitraten"          
[4] "ferf"               "co2C:nitraten"      "co2C:ferf"         
[7] "nitraten:ferf"      "co2C:nitraten:ferf"

lrt <- glmLRT(fit)
summary(decideTests(lrt))

       co2C:nitraten:ferf
Down                  986
NotSig              18833
Up                   1231

plotMD(lrt)

pval = 0.01
# et <- exactTest(y)
best = topTags(lrt, n = 10000)
best = best[best$table$FDR < pval, ]
# best = best[abs(best$table$logFC) > 1,]

# genes <- rownames(cpm(y, normalized.lib.sizes=TRUE))
DEgenes = rownames(best$table)
plotBCV(y)


abline(h = c(-1, 1), col = "blue")
abline(h = c(0), col = "red")

print(summary(decideTests(lrt, p.value = pval)))

       co2C:nitraten:ferf
Down                  556
NotSig              19705
Up                    789

# plots
plotMDS(y, top = 100)

Clustering sur les gènes DE

On utilise des modèles de mélange pour regrouper les gènes ayant des profils d’expression similaires dans les différentes conditions.

dataC <- data[DEgenes, ]
rowSums(dataC)

AT4G13180 AT4G15550 AT1G76680 AT2G39040 AT3G25250 AT1G29100 AT1G05660 AT1G78340 
    27625      6472     21867     18629      8410     15636      6423     19954 
AT2G29420 AT4G34131 AT1G17170 AT3G25730 AT1G12740 AT2G29490 AT4G12520 AT5G11920 
    14303     12224     27569     28345     16087     25883     21121     61112 
AT4G30420 AT4G22460 AT1G17860 AT4G15760 AT5G59540 AT1G17180 AT1G53830 AT2G29480 
    18662      9066     59045     12821     42467     25173     34734      3653 
AT1G13590 AT1G51830 AT1G79410 AT4G01870 AT5G26690 AT4G33790 AT4G19750 AT1G55920 
    17469     35374      4738     18652     12445     12613     31528     58830 
AT1G76690 AT1G12080 AT5G22890 AT4G19680 AT3G23560 AT3G11370 AT2G05380 AT4G34135 
    17508    102672     11373      7643     18332     13695       519     11128 
AT2G25810 AT2G47540 AT5G08790 AT5G63790 AT4G04460 AT1G78380 AT3G28345 AT4G19880 
    33515     15120     49996     39981     12163    170893     24028    139072 
AT3G09270 AT1G16370 AT3G20370 AT4G22070 AT4G14390 AT5G57540 AT4G31470 AT5G44130 
    44255      6917     49254      5244     12496      9931      7120     16659 
AT1G06120 AT2G41380 AT1G65680 AT2G43870 AT3G59140 AT2G33710 AT4G12510 AT1G65570 
    10592      1693     14040      3184      4259      2595      4003      8104 
AT3G09260 AT1G55790 AT4G04710 AT2G46660 AT1G13670 AT4G27860 AT1G54000 AT3G50480 
  1771228      5268      5121      4444     22163     32834    250063     29591 
AT2G25450 AT1G77450 AT3G18200 
   200340      2485     23513 
 [ reached getOption("max.print") -- omitted 1270 entries ]

conds <- str_split_fixed(colnames(dataC), "_", 2)[, 1]
run_pois <- coseq(dataC, conds = conds, K = 8:12, model = "Poisson", iter = 5, transformation = "none")

****************************************
coseq analysis: Poisson approach & none transformation
K = 8 to 12 
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1003.79063681148"
[1] "Log-like diff: 491.811396973711"
[1] "Log-like diff: 308.870088237903"
[1] "Log-like diff: 301.083062842708"
[1] "Log-like diff: 915.959558418981"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.00855398838491439"
[1] "Log-like diff: 0.00414549792929364"
[1] "Log-like diff: 0.00195388595650137"
[1] "Log-like diff: 0.000907633485260462"
[1] "Log-like diff: 0.000418364912153635"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.000245046326192266"
[1] "Log-like diff: 1.12953243309732e-05"
[1] "Log-like diff: 1.35195916683983e-07"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 74.2306168456969"
[1] "Log-like diff: 72.7178473602865"
[1] "Log-like diff: 6.7185891795049"
[1] "Log-like diff: 37.6742999938921"
[1] "Log-like diff: 99.4117559151978"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 112.146776250047"
[1] "Log-like diff: 141.097904817147"
[1] "Log-like diff: 105.378536868339"
[1] "Log-like diff: 138.9303123944"
[1] "Log-like diff: 177.562648327139"

coseq::plot(run_pois)

$logLike

$ICL

$profiles

$boxplots

$probapost_boxplots

$probapost_barplots

$probapost_histogram

summary(run_pois)

*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 12
ICL of selected model: -1647008
*************************************************
Number of clusters = 12
ICL = -1647008
*************************************************
Cluster sizes:
 Cluster 1  Cluster 2  Cluster 3  Cluster 4  Cluster 5  Cluster 6  Cluster 7 
        88        227        176        169        156         47         66 
 Cluster 8  Cluster 9 Cluster 10 Cluster 11 Cluster 12 
        39         43        183         56         95 

Number of observations with MAP > 0.90 (% of total):
1332 (99.03%)

Number of observations with MAP > 0.90 per cluster (% of total per cluster):
 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
 86        226       175       169       155       47        64       
 (97.73%)  (99.56%)  (99.43%)  (100%)    (99.36%)  (100%)    (96.97%) 
 Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
 39        42        178        56         95        
 (100%)    (97.67%)  (97.27%)   (100%)     (100%)    

clusters_per_genes <- clusters(run_pois)

ACP

On représente le clustering dans la plan principal d’une ACP

suppressMessages(library(ade4, warn.conflicts = F, quietly = T))
suppressMessages(library(adegraphics, warn.conflicts = F, quietly = T))

acp <- dudi.pca(log(dataC[, colnames(dataC) != "cluster"] + 0.1), center = TRUE, 
    scale = TRUE, scannf = FALSE, nf = 4)
summary(acp)

Class: pca dudi
Call: dudi.pca(df = log(dataC[, colnames(dataC) != "cluster"] + 0.1), 
    center = TRUE, scale = TRUE, scannf = FALSE, nf = 4)

Total inertia: 24

Eigenvalues:
     Ax1      Ax2      Ax3      Ax4      Ax5 
19.41941  3.31220  0.51028  0.23542  0.09744 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4     Ax5 
80.9142 13.8008  2.1262  0.9809  0.4060 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
  80.91   94.72   96.84   97.82   98.23 

(Only 5 dimensions (out of 24) are shown)

library("RColorBrewer")

dataC$cluster = clusters_per_genes[as.vector(rownames(dataC))]
s.corcircle(acp$co, xax = 1, yax = 2, fullcircle = FALSE, pback.col = "lightgrey")

adegraphics::s.class(acp$li, xax = 1, yax = 2, as.factor(dataC$cluster), labels = as.character(levels(as.factor(dataC$cluster))), 
    col = brewer.pal(n = 10, name = "Paired"), chullSize = 1, ellipseSize = 0, plabels.cex = 0.7, 
    pbackground.col = "grey85", main = "Clusters dans le plan principal", ylim = c(-9, 
        9))

s.corcircle(acp$co, xax = 2, yax = 3, fullcircle = FALSE, pback.col = "lightgrey")

adegraphics::s.class(acp$li, xax = 2, yax = 3, as.factor(dataC$cluster), labels = as.character(levels(as.factor(dataC$cluster))), 
    col = brewer.pal(n = 10, name = "Paired"), chullSize = 1, ellipseSize = 0, plabels.cex = 0.7, 
    pbackground.col = "grey85", main = "Clusters dans le plan principal", ylim = c(-9, 
        9))

s.corcircle(acp$co, xax = 4, yax = 2, fullcircle = FALSE, pback.col = "lightgrey")

adegraphics::s.class(acp$li, xax = 4, yax = 2, as.factor(dataC$cluster), labels = as.character(levels(as.factor(dataC$cluster))), 
    col = brewer.pal(n = 10, name = "Paired"), chullSize = 1, ellipseSize = 0, plabels.cex = 0.7, 
    pbackground.col = "grey85", main = "Clusters dans le plan principal", ylim = c(-9, 
        9))

s.corcircle(acp$co, xax = 4, yax = 3, fullcircle = FALSE, pback.col = "lightgrey")

adegraphics::s.class(acp$li, xax = 4, yax = 3, as.factor(dataC$cluster), labels = as.character(levels(as.factor(dataC$cluster))), 
    col = brewer.pal(n = 10, name = "Paired"), chullSize = 1, ellipseSize = 0, plabels.cex = 0.7, 
    pbackground.col = "grey85", main = "Clusters dans le plan principal", ylim = c(-9, 
        9))

Il semblerait que l’ACP détecte dans un premier temps (avec le premier vecteur principal) la valeur moyenne d’expression, puis l’expression diff induite par le fer.

On peut faire des représentations dans le plan des seconds et troisième axes principaux (le premier traduit le fer, le second la carence en nitrate).

Sur le cercle de corrélations dans le plan 2 et 3, on voit bien que lorsque la carence fer est là, on peut différencier un effet nitrate. Le CO2 est, lui encore difficile à identifier, le 4eme axe de l’ACP ne renseignant pas beaucoup…

Fonctions d’analyse pairwise

dualDE <- function(labels, pval = 0.01) {
    # data
    d <- data[, grepl(labels[1], colnames(data)) | grepl(labels[2], colnames(data))]
    y <- DGEList(counts = d, group = str_split_fixed(colnames(d), "_", 2)[, 1])
    y$samples
    
    # filtre
    keep <- filterByExpr(y)
    y <- y[keep, , keep.lib.sizes = FALSE]
    
    # normalisation
    y <- calcNormFactors(y)
    y$samples
    norm <- cpm(y, normalized.lib.sizes = TRUE)
    not_norm <- cpm(y, normalized.lib.sizes = FALSE)
    rd_genes = sample(rownames(norm), 500)
    heatmap(not_norm[rd_genes, ], main = "Sans normalisation")
    heatmap(norm[rd_genes, ], main = "Avec normalisation")
    
    # estimation of dispersion and tests
    y <- edgeR::estimateDisp(y)
    et <- exactTest(y)
    best = topTags(et, n = 10000)
    best = best[best$table$FDR < pval, ]
    # best = best[abs(best$table$logFC) > 1,]
    
    # genes <- rownames(cpm(y, normalized.lib.sizes=TRUE))
    DEgenes = rownames(best$table)
    plotBCV(y)
    
    plotMD(et, p.value = pval)
    abline(h = c(-1, 1), col = "blue")
    abline(h = c(0), col = "red")
    print(summary(decideTests(et, p.value = pval)))
    
    # plots
    plotMDS(y, top = 100, col = rep(1:2, each = 3))
    logcpm <- cpm(y, log = TRUE)
    heatmap(logcpm[DEgenes, ])
    print("ON SELECTIONNE AU TOTAL :")
    print(paste(length(DEgenes), "genes with fdr < ", pval))
    
    results <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", 
        "description", "external_gene_name", "entrezgene_id"), values = DEgenes, 
        mart = mart)
    results <- results[!rownames(results) %in% which(duplicated(results$ensembl_gene_id)), 
        ]
    
    rownames(results) = results$ensembl_gene_id
    r = results[DEgenes, ]
    r$pval = best$table$PValue
    r = r[order(r$pval), ]
    return(r)
}


OntologyProfile <- function(r) {
    ego <- enrichGO(gene = r$entrezgene_id, OrgDb = org.At.tair.db, ont = "BP", pAdjustMethod = "BH", 
        pvalueCutoff = 0.01, qvalueCutoff = 0.05, readable = TRUE)
    simpOnt <- simplify(ego, cutoff = 0.7, by = "p.adjust", select_fun = min)
    simpOnt@result$Description
    
    print(barplot(simpOnt, showCategory = 40, font.size = 5))
    print(dotplot(simpOnt, showCategory = 40, font.size = 5))
    print(emapplot(simpOnt, layout = "kk"))
}



mat <- matrix(nrow = length(unique(conds)), ncol = length(unique(conds)))

comparables = c()
for (c1 in unique(conds)) {
    for (c2 in unique(conds)) {
        nb_diffs = 0
        for (i in seq(1, 3)) {
            if (str_split(c1, "")[[1]][i] != str_split(c2, "")[[1]][i]) {
                nb_diffs = nb_diffs + 1
            }
        }
        if (nb_diffs == 1) {
            comparables = c(comparables, paste0(c1, "-", c2))
        }
    }
}

getGEgenesNumber <- function(labelsConc, pval = 0.01) {
    # data
    labels = c(str_split_fixed(labelsConc, "-", 2)[, 1], str_split_fixed(labelsConc, 
        "-", 2)[, 2])
    print(labels)
    d <- data[, grepl(labels[1], colnames(data)) | grepl(labels[2], colnames(data))]
    y <- DGEList(counts = d, group = str_split_fixed(colnames(d), "_", 2)[, 1])
    y$samples
    
    # filtre
    keep <- filterByExpr(y)
    y <- y[keep, , keep.lib.sizes = FALSE]
    
    # normalisation
    y <- calcNormFactors(y)
    y$samples
    
    # estimation of dispersion and tests
    y <- edgeR::estimateDisp(y)
    et <- exactTest(y)
    best = topTags(et, n = 10000)
    best = best[best$table$FDR < pval, ]
    # best = best[abs(best$table$logFC) > 1,]
    
    # genes <- rownames(cpm(y, normalized.lib.sizes=TRUE))
    DEgenes = rownames(best$table)
    return(length(DEgenes))
}

Effet Nitrate

library(data.table)


labels <- c("cnF", "cNF")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       cNF-cnF
Down       758
NotSig   17187
Up        2054

[1] "ON SELECTIONNE AU TOTAL :"
[1] "2812 genes with fdr <  0.01"

# kable(r)

kable(r[r$description %like% "nitra", ])

kable(r[r$external_gene_name %like% "NRT", ])

kable(r[r$external_gene_name %like% "NR", ])

# NRT2.1
getExpression("AT1G08090")

[1] "all"

getExpression("AT1G08090", conds = c("cnF", "cNF"))

[1] "cnF" "cNF"

# NAR2
getExpression("AT5G50200")

[1] "all"

getExpression("AT5G50200", conds = c("cnF", "cNF"))

[1] "cnF" "cNF"

"AT5G50200" %in% r$ensembl_gene_id

[1] FALSE

OntologyProfile(r)

labels <- c("CnF", "CNF")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       CNF-CnF
Down       512
NotSig   18162
Up        1259

[1] "ON SELECTIONNE AU TOTAL :"
[1] "1771 genes with fdr <  0.01"

# kable(r)

kable(r[r$description %like% "nitra", ])

getExpression("AT1G08090", conds = labels)

[1] "CnF" "CNF"

getExpression("AT5G50200", conds = labels)

[1] "CnF" "CNF"

OntologyProfile(r)

Effet CO2

labels <- c("cNF", "CNF")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       CNF-cNF
Down        46
NotSig   19757
Up          18

[1] "ON SELECTIONNE AU TOTAL :"
[1] "64 genes with fdr <  0.01"

kable(r)

kable(r[r$description %like% "oxid", ])

OntologyProfile(r)

labels <- c("cnF", "CnF")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       CnF-cnF
Down       666
NotSig   18680
Up         746

[1] "ON SELECTIONNE AU TOTAL :"
[1] "1412 genes with fdr <  0.01"

OntologyProfile(r)

labels <- c("cNf", "CNf")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       CNf-cNf
Down      1040
NotSig   18086
Up        1232

[1] "ON SELECTIONNE AU TOTAL :"
[1] "2272 genes with fdr <  0.01"

OntologyProfile(r)

labels <- c("cnf", "Cnf")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       Cnf-cnf
Down       603
NotSig   18788
Up        1031

[1] "ON SELECTIONNE AU TOTAL :"
[1] "1634 genes with fdr <  0.01"

OntologyProfile(r)

Effet Fer

labels <- c("cNf", "cNF")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       cNF-cNf
Down      3497
NotSig   12986
Up        3923

[1] "ON SELECTIONNE AU TOTAL :"
[1] "7420 genes with fdr <  0.01"

OntologyProfile(r)

labels <- c("cNf", "CNf")
r <- dualDE(labels, pval = 0.01)

Design matrix not provided. Switch to the classic mode.

       CNf-cNf
Down      1040
NotSig   18086
Up        1232

[1] "ON SELECTIONNE AU TOTAL :"
[1] "2272 genes with fdr <  0.01"

OntologyProfile(r)

En résumé : meeting Antoine Sophie

nbDe <- sapply(comparables, getGEgenesNumber)

[1] "cNF" "cnF"
Design matrix not provided. Switch to the classic mode.
[1] "cNF" "CNF"
Design matrix not provided. Switch to the classic mode.
[1] "cNF" "cNf"
Design matrix not provided. Switch to the classic mode.
[1] "cnF" "cNF"
Design matrix not provided. Switch to the classic mode.
[1] "cnF" "CnF"
Design matrix not provided. Switch to the classic mode.
[1] "cnF" "cnf"
Design matrix not provided. Switch to the classic mode.
[1] "CNF" "cNF"
Design matrix not provided. Switch to the classic mode.
[1] "CNF" "CnF"
Design matrix not provided. Switch to the classic mode.
[1] "CNF" "CNf"
Design matrix not provided. Switch to the classic mode.
[1] "CnF" "cnF"
Design matrix not provided. Switch to the classic mode.
[1] "CnF" "CNF"
Design matrix not provided. Switch to the classic mode.
[1] "CnF" "Cnf"
Design matrix not provided. Switch to the classic mode.
[1] "cNf" "cNF"
Design matrix not provided. Switch to the classic mode.
[1] "cNf" "cnf"
Design matrix not provided. Switch to the classic mode.
[1] "cNf" "CNf"
Design matrix not provided. Switch to the classic mode.
[1] "cnf" "cnF"
Design matrix not provided. Switch to the classic mode.
[1] "cnf" "cNf"
Design matrix not provided. Switch to the classic mode.
[1] "cnf" "Cnf"
Design matrix not provided. Switch to the classic mode.
[1] "Cnf" "CnF"
Design matrix not provided. Switch to the classic mode.
[1] "Cnf" "cnf"
Design matrix not provided. Switch to the classic mode.
[1] "Cnf" "CNf"
Design matrix not provided. Switch to the classic mode.
[1] "CNf" "CNF"
Design matrix not provided. Switch to the classic mode.
[1] "CNf" "cNf"
Design matrix not provided. Switch to the classic mode.
[1] "CNf" "Cnf"
Design matrix not provided. Switch to the classic mode.

df <- data.frame(Comparisons = comparables, DEGNumber = nbDe)
df <- df[rownames(df)[which(!duplicated(df$DEGNumber))], ]
kable(df[order(-df$DEGNumber), ])

library(reshape2)

p <- ggplot(data = df, aes(x = Comparisons, y = DEGNumber, fill = Comparisons)) + 
    geom_dotplot(binaxis = "y", stackdir = "center") + theme(axis.text.x = element_text(angle = 320, 
    hjust = 0, colour = "grey50"), plot.title = element_text(size = 14, face = "bold")) + 
    ggtitle(paste("Number of DE genes for every condition comparisons"))
p